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Abstract 

Dyson-Schwinger equations are important tools for non-perturbative analyses of 
quantum field theories. For example, they are very useful for investigations in quan- 
tum chromodynamics and related theories. However, sometimes progress is impeded 
by the complexity of the equations. Thus automatizing parts of the calculations will 
certainly be helpful in future investigations. In this article we present a framework 
for such an automatization based on a C++ code that can deal with a large number 
of Green functions. Since also the creation of the expressions for the integrals of 
the Dyson-Schwinger equations needs to be automatized, we defer this task to a 
Mathematica notebook. We illustrate the complete workflow with an example from 
Yang-Mills theory coupled to a fundamental scalar field that has been investigated 
recently. As a second example we calculate the propagators of pure Yang-Mills 
theory. Our code can serve as a basis for many further investigations where the 
equations are too complicated to tackle by hand. It also can easily be combined 
with DoFun, a program for the derivation of Dyson-Schwinger equationsj^ 

Keywords: Dyson-Schwinger equations, correlation functions, quantum field 
theory 



PROGRAM SUMMARY 

Program Title: CrasyDSE 
Version: 1.1.0 

Licensing provisions: CPC non-profit use license 
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PACS: 11.10.-z,03.70.+k,11.15.Tk 
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in C++. This code uses structures to handle large numbers of Green functions. 

Unusual features: Provides a tool to convert Mathematica expressions into C++ ex- 
pressions including conversion of function names. 

Running time: Depending on the complexity of the investigated system solving the 
equations numerically can take seconds on a desktop PC to hours on a cluster. 



1. Introduction 

Strongly coupled field theories play an essential role in the physical description 
of nature. Both established theories like quantum chromodynamics and conjectured 
ones like technicolor theories cannot be fully understood without non-perturbative 
methods. Typical approaches include Monte-Carlo simulations on a discretized 
space-time or functional equations. Functional renormalization group equations, 
see, e. g., [IHl], Dyson-Schwinger equations, see, e. g., [5H8], and the n-PI formal- 
ism, see, e. g., 0, belong to the second group. Their advantages are well appreciated 
and they provided many new insights. 

In this article we will focus on Dyson-Schwinger equations (DSEs) and propose 
a concrete way to handle them when they become too complex to be treated by 
hand alone. DSEs consist of a system of coupled integral equations which relate 
different Green functions. Since there are infinitely many Green functions there 
are also infinitely many DSEs. Unfortunately no subset of these equations forms a 
closed system so that we have to deal with an infinitely large system of equations. 
Naturally one hopes that only a (small) finite number of Green functions is relevant 
and looks for truncations capturing the most important features of the theory. Of 
course, in order to check the validity of such an approach one should test the influ- 
ence of neglected Green functions. However, this is often very tedious. On the other 
hand there are also theories where it is known that current truncation and approx- 
imation schemes and available methods are insufficient and need to be extended. 
For example, standard truncations restrict the DSEs to one-loop diagrams [5l [H [TOl - 
\n\ but Yang-Mills theories in the maximally Abelian gauge require the inclusion 
of two-loop diagrams in order to be consistent in the non-perturbative regime jl3] . 
Consequently the present technical methods have to be improved. 

The reason why more sophisticated truncations or more complicated theories 
require so much more effort is mainly that the length and complexity of the explicit 
expressions increase considerably with the number of interactions and the number 
of external legs. Also the numbers of dressing functions and diagrams grow for 
higher Green functions. We will illustrate this below explicitly with the example of 
Yang-Mills theory coupled to a scalar: We will see that extending a simple truncation 
beyond the propagators by dynamically including the vertex between the gauge field 
and the scalar triples the number of dressing functions to be calculated and requires 
five times as many loop integrations. Furthermore, the corresponding integration 
kernels are substantially more complicated than the first one. Seeing such complexity 
arise from such a simple extension we felt it was time to think about automatizing 
this process. This seems even more necessary since computing time is no longer as 
restrictive as it was ten years ago. For example, fourteen years ago the first solution 
of the DSE system of Yang- Mills theory that was complete at the propagator level 
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[TOl [HI [15] relied on an angle approximation and took several hours. Nowadays it is 
possible to do it with the full momentum integration in a few minutes. Thus more 
complicated truncations and theories are definitely doable. However, right now one 
has to invest much time in deriving DSEs and implementing them. In a sense we 
fell behind the possibilities today's computers offer and we think we should try to 
change this and find means that allow us to focus more on the physical rather than 
the technical problems. 

The technical part of investigating a theory numerically with DSEs consists of 
two main steps: First, the equations have to be derived. Second, one has to imple- 
ment them in a numeric code. A tool that assists in the first part is already available 
with the Mathematica [12] application DoFun [T7]. Here we present a generic nu- 
meric code that can serve as a basis for the second step. CrasyDSE (Computation 
of Rather lArge SYstems of DSEs) is capable of dealing with a high number of 
Green functions and their dressing functions. Furthermore it provides several pre- 
defined integration routines and numerical approximation techniques. It can also 



be extended to multi-core environments (see comments in section 4.2) and finite 



temperature (see comments in section 4.1). Part of CrasyDSE is the Mathematica 



package CrasyDSE. m to generate C++ expressions for the kernels. This first of 
all alleviates the generation of the code tremendously and reduces human errors 
and secondly is the easiest way to transform the notation of the user into the nota- 
tion of CrasyDSE. Note that the functions of the package can deal with all regular 
Mathematica input and do not rely on DoFun. In order to use the package, the 
file CrasyDSE. m has to be copied from main_Mathematica to a place where Math- 
ematica can find it. We suggest to copy it to the subdirectory Applications of 
the Mathematica user directory (SUserBaseDirectory/ApplicationsJ^ Now the 
package can be loaded with <<CrasyDSE'. 

In the following we will describe the general procedure to solve DSEs in section [2j 
The numerical problem is formulated in section [3] and section |4] contains details 
on the provided routines to solve DSEs. Sees. [5] and [6] explain the application of 
CrasyDSE using as examples the calculation of two DSEs of a scalar field coupled to 
Yang-Mills theory and the calculation of the propagators of pure Yang-Mills theory. 
Finally, we give a summary and an outlook in section [7} In three appendices we 
provide details and summaries on functions and variables of CrasyDSE. 

2. Solving Dyson-Schwinger equations 

Our approach to solving DSEs can be separated into three parts as illustrated 
in Fig. [T| 

1. Derive the equations from the given action by the method of choice. If not 
done in Mathematica, enter them into Mathematica. 

2. Use the Mathematica package CrasyDSE. m to generate the C++ files with the 
kernels. Alternatively, in simple cases one can write the kernel files manually. 

3. Use the kernel files with the C++ code of CrasyDSE to solve the DSE numer- 
ically. 

We illustrate the last steps with two examples in sections |5] and [6j 



^On a Unix machine this is typicahy ~/ .Mathematica/ Applications, on Windows it is 
User\ Application Data\Mathematica\ Applications. 
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Figure 1: Schematic workflow. 



For the first step, the derivation of the DSEs, we recommend the Mathematica 
application DoFun (Derivation Of FUNctional equations) [TTj. Its predecessor is 
the Mathematica package DoDSE (Derivation Of DSEs) jSl [18] which was of great 
help in the investigation of big systems of DSEs like that of the maximally Abelian 
gauge [13] and the Gribov-Zwanziger action [20] . The calculation of some in- 
frared properties in the maximally Abelian gauge would even have been impossible 
without automatization due to the huge number of terms [TSl |2I] • Later on DoDSE 
was considerably extended and the derivation of functional renormalization group 
equation^ was included [17], thus it was renamed to DoFun. However, note that 
CrasyDSE works completely independent of DoFun. 

The second step consists in making the Mathematica expressions accessible for 
C++. We chose to write our own functions that generate complete C++ files. 
Thus we maintain as much control as possible over the process and make it more 
transparent for the user. However, in principle it would also be possible to let 
Mathematica and C++ interact in a more direct way via MathLink. All the necessary 
functions to create the C++ files are included in the package CrasyDSE. m, but the 
notebook from which it is created is also provided so that the user has direct access 
and can most easily adapt code if required. 

Finally, after all kernels have been written into C++ files, one uses the pro- 
vided C++ modules to solve the equations. Typical initial work includes defining 
model parameters, defining the required Green functions and their dressing func- 
tions, choosing integration routines and defining the renormalization procedures. 
The way dressing functions are defined is thereby quite arbitrary: One can use closed 
expressions, for example, from fits to lattice data, interpolations or expansions in 
sets of polynomials. In order to help with starting calculations with CrasyDSE we 
provide several examples with the main code. 

We want to stress that CrasyDSE can not be considered a black box. In order 
to successfully use it, the user has to understand many of the employed routines 
and adapt them if required. CrasyDSE is merely a framework for solving DSEs 
that provides structures to handle Green functions and their DSEs and modules 
to perform the most basic steps like integration. However, as every problem has 
its own intricacies the user still has to implement many specific functions, e. g., 
extrapolation functions for dressing functions. 



^Recently a similar program to CrasyDSE has become available for functional renormalization 
group equations with the program FlowPy [55] . 
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3. General formulation of the numerical problem 



As already mentioned DSEs form an infinitely large set of equations. For nu- 
meric calculations we take a subset of these equations, but they will always depend 
on Green functions whose DSEs are not part of the truncation. Depending on the 
details of our truncation scheme we can either provide expressions for them as ex- 
ternal input or drop diagrams containing such Green functions. Furthermore, every 
Green function can consist of several dressing functions and before we can do any 
calculation we have to project the DSEs such that we deal with scalar integrals. 
Sometimes it is not possible or not feasible to project directly onto the dressing 
functions and additionally a linear system of equations has to be solved to get re- 
sults for them. But let us for now assume for simplicity that it is possible to project 
directly onto the dressing functions. 

Every Green function depends on several momenta. The dressing functions, 
however, depend only on a reduced number of variables. For example, a two-point 
function depends on two external momenta. Momentum conservation reduces these 
to one momentum. The dressing function (s) of such a Green function depend only 
on the remaining momentum squared, i. e., one variable instead of four. For a 
three-point function there are two independent momenta and the dressing functions 
depend on three variables. In a slight abuse of language we call the variables of the 
dressing functions external momenta. The number of variables is denoted as the 
dimension of external momenta. Later we will also encounter internal momenta. 
These are the remaining loop momentum variables after trivial integrations have 
been performed. 

In the following dressing functions are denoted by where g denotes the 

Green function to which it belongs and i labels the dressing functions of a Green 
function. Xg G Qg represents the external momenta. We have to solve the following 
integral equations: 

/ 

A^^' : fig C M"'^ ^ M , 

where ^^are ^'^^ bare dressing functions (if non-zero and possibly including a 
renormalization constant). The sum over / denotes contributions from different 
graphs and Z^'*'^ are the renormalization constants of the bare n-point function 
of a given graph. The integration is over the loop momenta y. The dimensions 
of external and internal momenta are indicated by dg and di, respectively. The 
F/ (?/, Xg, {A} , {Amodoi}) denote the kernels of the integrals. They depend on the 
internal and external momenta explicitly and via several Green functions also im- 
plicitly. Some of them, the {A}, are a dynamic part of the truncation, whereas 
others, the Amodei, are given by external input. 

Before solving this system of equations numerically the following steps are re- 
quired: 

1. The dressing functions A^'* have to be approximated, e. g., by discretization of 



the argument or expansion in an orthogonal set of functions (see section 4.1). 



2. If the integrals are divergent, one needs a regularization prescription, e. g., a 



sharp cutoff in \y\ or BPHZ [231425] (see section 5.1). 
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3. Expressions for the dressings Amodei need to be provided (see sections 4.3 and 



5.1) 



4. A renormalization procedure needs to be defined to fix the renormahzation 
constants Z^'^'^ (see sections 4.3 and 5.1 and Appendix B). 



After this is settled, the integrals can be evaluated with the provided quadratures 
(see section 4.2). A solution can be found, for example, by fixed point iteration or 



Newton's method, see, e. g., [261428] . Both these solution methods are implemented 



in CrasyDSE, see section 4.3 



4. Implementation in C+ + 

Solutions to the different problems stated in section [3] are implemented in C++ 
in three modules: 

• dressing, cpp/hpp, 

• quadrature, cpp/hpp, 

• DSE. cpp/hpp, 

where dressing, cpp/hpp uses the structure dse from DSE. cpp/hpp. 

Additionally some simple general functions are stored in function. cpp/hpp. All 
these files can be found in the directory main. The provided examples are located 
in separate directories in examples. They are initialized and called in the files 
sphere.main.cpp, interp_main.cpp, scalar_main.cpp and YM4 domain, cpp. A very 
basic Makefile for Unix using the g++ GNU compiler is provided with each of the 
examples. 

The provided modules are as self-contained as possible and can be arbitrar- 
ily extended, e. g., by including further interpolations in dressing, cpp/hpp, adap- 
tive integration algorithms in quadrature. cpp/hpp or additional solving strategies in 
DSE. cpp/hpp. 

4.1. Approximation (dressing. cpp/hpp) 

All functions and parameters relevant for approximating the dressing functions 
are referenced to or stored in the structure struct dse defined in DSE.hpp. A 



summary together with the corresponding objects in a DSE is given in table C.6 

The expansion coefficients for the int dim_A dressing functions are contained in 
the array double *A. We have int dim_x external variables, and the numbers of ex- 
pansion coefficients for each of them are saved in the array int *n_A. In the case of 
linear interpolation the interpolation points have to be stored in the array double 

diin_x-l 

*x_A. The total number of expansion coefficients is int ntot_A:= Yl ii_A[i]. 

i=0 

Since future extensions for non-zero temperature calculations require also discrete 
variables, i. e., Matsubara frequencies, part of the variables can be considered as 
integer numbers z of dimension int dim_mat. Currently non-zero temperature is 
not fully implemented and thus dimjnat should always be set to zero. The al- 
location/deallocation of the arrays A and x_A is done with void init_A_xA(void 
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*dse_param) /void dealloc_A_xA(void *dse_param)|^ A minimal dse structure 
can be initialized with void init_def ault_dse(void *dse_parain) . More detailed 



information on A and x_A can be found in [Appendix A 



In dressing, cpp/hpp two ways to express the dressing functions are provided: 

• linear interpolation (Lin_gen_dress_interp), 

• expansion in Chebyshev polynomials (Cheb_gen_dress_interp). 

In the case of a linear interpolation A contains the function values on a rectilinear 
grid stored (together with the discrete arguments) in x_A whereas in the case of an 
expansion in Chebyshev polynomials the expansion coefficients are stored in A and 
only the discrete arguments need to be stored in x_^[^ For a Chebyshev expansion, 
as introduced for dressing functions of Green functions in [26], it is possible to 
either transform the standard interval [—1,1] linearly or logarithmically to the actual 
domain of interpolation via setting int cheb_traf o [i] to 1 or 2, respectively, where 
i denotes the external variable. By default it is set to 1 in init_A_xA. Furthermore, 
with int cheb_f unc_traf o [i] one can expand the logarithm of a function |26] 
instead of the function itself by changing its default value 0, set in init_A_xA, to 1. 



i denotes here the dressing function. For an example see section 6.2 

The provided interpolation algorithms work only as long as the arguments x of 
the dressing function are inside the user-defined domain Qg. To determine if x E 
Qg the user-defined function void def _domain(double *x, void *dse_param) is 



called. For details of how to construct this function see Appendix B If x ^ ilg, the 
user has to provide a function for extrapolation via double interp_off domain (int 
*pos, double *x, int iA, void *dse_param). Details can again be found in 



Appendix B but the gist is that the array pos knows on which side of the allowed 
interval the external variable Xi lies: If Xi < ai or Xi > bi, pos[i] is 1 or 2, 

respectively. 

The correct initialization and definition of the required parameters and functions 
is illustrated by the example interp_main.cpp interpolating three functions linearly 
and with Chebyshev polynomials. These functions have two discrete and three 
continuous variables where the domains of the continuous variables depend on the 
values of the discrete variables. 

4-2. Integration (quadrature. cpp/hpp) 

All functions and parameters relevant for integration are referenced to or stored in 
the structure struct quad. It is defined in quadrature. hpp and allocation and deallo- 
cation is done with void init_quad(void *quad_param) and void dealloc_quad(void 
*quad_param) , respectively. An overview of the members of quad relevant to the 



user and their usage in the context of DSEs is given in table C.7 

This module provides the means to integrate nint integrals of dimension dim. 
Any integral is split into three parts: a constant factor independent of any vari- 
ables, a Jacobian and the remaining integrand. These three parts have to be given 



^In order not to overload the text we refrain in most cases from a detailed explanation of all 
arguments and only give them for reference. Details on the meaning of the arguments can be found 
directly in the code where all of them are explained in the function descriptions. 

^Note that when using the DSE solving routines also for Chebyshev interpolation the continuous 
external momenta have to be stored in x_A using the function Cheb_init_coiit_xA described in 

Appendix "xn 



7 



variable 1 with 

n_part[0]=3 integration regions: 



variable 2 with 

n_part[l]=2 integration regions: 



nodes_part[0] 

type[0] 

traf[0] 



nodes_part[l] 

type[l] 

traf[l] 



nodes_part[2] 

type[2] 

traf[2] 



nodes_part[3] 

type[3] 

traf[3] 



nodes_part[4] 

type[4] 

traf[4] 



Figure 2: Example of how integration regions and related variables are used: For two integration 
variables the array n_part contains the information how many integration regions there are for 
each variable. The resulting integration regions are then numbered consecutively. The number of 
integration points, the integration type and the transformation type of each region are stored in 
the arrays nodes_part, type and traf . 



by the three functions void coeff (double *coef f icients , void *int_param) , 
double Jacob (double *x) and void integrand (double *erg, double *x, void 
*int_param). They have to be defined by the user, but integrand and coeff are 
usually created with the CrasyDSE Mathematica notebook. Further details on these 



functions can be found in Appendix B 



The boundaries of the integrals are defined in the function void boundary (double 
*bound, double *x, int idim, void *int_param). Details on its required con- 



tents are provided in Appendix B For now it suffices to say that it defines the 
domains of the integration variables yi, where the inner integration bound- 

aries may depend on the outer integration variables: 

bg bi{yo) bdim-iCj/o,-- -11/(1111.-2) 

dyo j dyi--- j (i?/dim-i- (2) 

Boundaries can also depend on the external momenta. However, in this case the 
integration routine becomes slower. If this is required one has to set int bound_type 
to 1. The default value, set in init_quad, is 0, i. e., the boundaries must not 
depend on the external momenta. As an example where the boundaries depend on 
the external momenta one can consider the integration of the Yang- Mills system in 
section [6l 

The integration of every variable yi can be split into several parts which may 
increase the precision of the results, see, e. g., The number of integration regions 
for the i-th integration variable is given by int n_part [i] . Each region is labeled 
here by j and the number of its integration points is set by int nodes.part [j] . 
Fig. [2] illustrates the numbering of integration regions. For each region a quadrature 
rule has to be chosen via setting int type [j] to a number corresponding to one of 
the quadrature rules given in table [T| where also necessary parameters are indicated. 
The quadrature rules are defined on the interval [—1, 1], which can be transformed 
with various functions to the actual integration interval defined in boundary. This 
works by setting int traf [ j ] to one of the values indicated in table |2j 

Finally we want to draw attention to the fact that the integration is the most 
costly part of solving DSEs. Therefore a parallelization of the program is most 
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J -J- 1 

quadrature rule 


typeLj] 


parameters 


Gauss-Legendre 





none 


Gauss-Chebyshev type two 


1 


none 


Fejers second rule 


2 


none 


Double exponential 


3 


int param [0] : stepsize 



Table 1: Currently implemented quadrature rules, the corresponding values of type[j] and the 
required parameters. Details can be found in quadrature, cpp. 



transformation rule 


traf [j] 


none 





linear 


1 


logarithmic 


2 


modified logarithmic 


3 


modified logarithmic [2Sj 


4 



Table 2: Currently implemented transformation rules from [—1,1] to the actual integration interval 
and the corresponding values for traf [ j ] . Details on these transformations can be found in the 
function nw_trafo of quadrature. cpp. 



efficient in the function void integrate (double *erg, void *quad_param, void 
*int_param) of the quadrature module, where the loop over the internal or external 
momenta can be distributed to several cores. This is not implemented but a user 
familiar with parallelization should be able to extend the program in this direction 
without too much effort. 

To summarize the user has to provide the following functions: 

• integrand: Defines the integrand of the integral. Usually this will be the 
kernel function created by the Mathematica notebook. 

• coef f : Constant factor. Usually it is created by the Mathematica notebook. 

• Jacob: Defines the Jacobian of the integral measure. 

• boundary: Initializes the boundaries Oq, b^, adim-i(l/0) • • • ; Z/dim-2) and 
bdim-i{yo, ■ ■ ■ ,ydim-2), where each boundary can depend on previous integra- 
tion variables. 

The correct initialization and definition of the needed parameters and functions 
is illustrated by a simple example. In sphere^main.cpp the function 

sphere_integrand : , (3) 

times the Jacobian r is integrated over 

R 2n s/W^ 

dz I d<\) [ dr. (4) 



'R 
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Furthermore, we exemplify here the use of external parameters. For solving DSEs 
the external parameters are normally the external momenta. However, the inte- 
gration routine can handle also other cases of external parameters. In general one 
can define int nint_para different parameters initialized in void init_para(int 
i , void *int_param) for which the integral is performed by calling the integration 
routine once. For the standard application in DSEs these functions are set auto- 
matically as required by the function init_A_xA. In the example sphere ^main.cpp 
another possibility is demonstrated. Here init_para is set to sphere_init_para 
which sets the external parameters as multiplicative factors for the integrals. 

4-3. Solving DSEs (DSE.cpp/hpp) 

Assuming that together with the quadrature a proper regularization has been 
chosen all the integrals in Eq. ([T]) are known and we are left with the task of solving 
the given integral equations including a proper renormalization. 

Every Green function present in our truncated set of DSEs is represented by 
its own structure dse which contains all the necessary functions and parameters in 
order to evaluate its dressing functions via the function dress. This is already all 
one needs to initialize the modeled Green functions, whereas those we are going to 
solve for need additional information in their structures. A summary of all variables 



and functions relevant for the user is given in table C.8 The fact that the equa- 
tions are coupled and the iteration of one dressing function needs information about 
dressing functions of int n.otherGF other Green functions is handled by struct 
dse *otherGF which contains copies of the needed structures, called when evalu- 
ating the integration kernels. Therefore it is necessary to allocate the arrays in all 
other Green functions before copying them to otherGF such that the correct pointer 
addresses are available. Only variables which are not supposed to be changed during 
the iteration procedure will be copied by value. Additionally some model param- 
eters might be needed by the (modeled) dressing functions which are pointed to 
by every structure dse via struct mod, defined by the user. For an example see 
Fig. [3| where also the quad structures are indicated which contain specifics on the 
integration routines needed to evaluate the loop integrals in the graphs. 

Focusing now on one specific Green function these graphs are grouped into int 
n_looporder contributions where all int n_loop [i] members of one of the groups 
can be evaluated using the same quadrature struct quad Q [i] . The evaluation 
of the self-energy for the int dim_A dressing functions and int ntot_A different 
array points Xg is then performed by calling void self energy (void *dse_param) . 
The result is stored in double *self _A. It will be used in the user-defined function 



void renormCvoid *DSE) to obtain the new dressing functions, see Appendix B 
for details. 

4-3.1. Iteration 

We implemented a fixed point iteration solution technique, i. e., the system is 
solved by calculating 
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three-gluon 
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quadrature 



scalar 
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scalar 
propagator 



quadrature 



Figure 3: Schematic illustration of structures defined in the code and their dependencies for the 
example of Yang-Mills theory coupled to a scalar field considered in section [5j 



from the previous dressing functions ^^fc)- The calculation of the new dressing 

functions ^[^+1) is performed in several steps, where it might be necessary that a 
subset of Green functions is iterated till convergence for every "meta-iteration" of 
the full set of equations. 

Before starting the iterations the initial dressings A^q^ of every Green function 
have to be set via void init_dress (void *dse_param) . The last step in every iter- 
ation is then to renormalize such that the ^('^+1) have the correct values /derivatives 
at the renormalization scale (s) /j, via int renorm_n_parain parameters which may 
need some initialization in void init_renormparain(void *dse_param). In the ex- 
ample of section |5] renormalization will be done by fixing the int renorm_n_Z renor- 
malization constants Z^^^-^'' accordingly in the function void renorm(void *DSE). 

Note that renormalization constants also appear in the bare Green functions A^^*^ 
but one could employ subtracted equations to drop them. 

The iteration of a single Green function is performed in void solve_iter (void 
*dse_param, int output) where different stopping criteria (absolute difference 
double epsabs, relative difference double epsrel and maximum number of itera- 
tions int maxiter) are available. Several Green functions can be united in an ar- 
ray which can be passed to void meta_solve_iter (struct dse *DSE, int ndse, 
double epsabsstop, double epsrelstop, int iterstop, int output). It has 
the same stopping criteria as solve_iter and in every meta-iteration solve_iter 
is called for every Green function. This allows different relative iterations, e. g., 
all Green functions are iterated once for every meta-iteration or a subset of Green 
functions is iterated till convergence while another subset is iterated only once. 

To get an idea of the efficiency of our code we compared the calculation from 
section [5] with an independently created code that was optimized for this problem 
[5U] . In general the difference in time depends on how much optimization is possible 
in the given problem. In the present case the time difference for one iteration was 
less than a factor of 2. 
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4-3.2. Newton's method 

Another method for solving DSEs is based on Newton's method to solve a non- 
linear system of equations. It was already used in many DSE calculations, see, for 
instance, [2SH2H1 EI]- For this method the system of DSEs is rewritten into the 
following form: 



(6) 



where i labels the dressing functions of all DSEs and Xk denotes the external mo- 
menta. We assume here that the dressing functions A^{xk) are expanded in a set of 
basis functions. The corresponding expansion coefficients are the unknown variables 
where j labels the polynomials. The goal is to find those values for c^*'-'-' that 
make all iJ^*'^) vanish. Newton's method yields new coefficients by the following 
formula: 



i',k 

where the Jacobian J is given by 

The backtracking parameter A can be used to optimize this step by choosing an 
appropriate value between zero and one. The determination of A can be subject of 
sophisticated algorithms, see, for example, [28]. Here, however, we simply cut A in 
half if the norm of the new is not smaller than that of the old one. If the 

starting functions are well chosen this is sufficient for the example of section |6j This 
procedure is repeated until the norm of the vector E'^^''^'* drops below a given value 
or a maximal number of iterations is reached. A single iteration step takes here 



much longer than for the fixed point iteration described in section 4.3.1 because the 
calculation of the Jacobian is rather expensive. In principle the derivatives required 
to get J can be done directly, but here Broyden's method is used which defines an 
approximate Jacobian by a simple forward differentiation with small step size h: 

'^{i',A:),approx ' 

where the notation E^"^' '^\c^^'^^ + h) means that is calculated with the coef- 

ficient c*^*'-'-' changed by /i, whereas E^'^ '''\c^^'^^) refers to the original E^'^ This 
prescription proved very reliable for the example treated in section [6] 

Newton's method is implemented in the function void solve_iter_secant (struct 
dse *DSE, int ndse, int maxiter, double eps_E, int output). As arguments 
it takes the array of dses DSE, the length of this array ndse, the maximal number of 
iterations maxiter, the stopping value for '^i^E^^''^^ eps_E and an integer number 
output which determines if intermediary output should be printed to the screen 
(1) or not (0). Of course Newton's method can be combined with the fixed point 
iteration. For example, one could solve two of three DSEs with Newton's method 
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and then iterate the third one with solve_iter. 



5. Solving the gap and vertex equations of Yang-Mills theory coupled to 
a scalar field 

In this section we describe how to solve a truncated set of DSEs of Yang-Mills 
theory coupled to a scalar field. We will first give a short overview of the employed 
truncation and then explain how to solve it with CrasyDSE. 

The derivation of the DSEs is not discussed here, but we provide details in the 
Mathematica notebook DoFun^YM+Scalar.nh. The results of this notebook form 
the basis on which we create the functions for the C++ code. The corresponding 
steps are contained in a second notebook, CrasyDSE_YM+scalar.nb. It describes 
how the expressions of the DSEs have to be modified so they can be used as input 
for CrasyDSE. In a second part all required definitions are provided and the kernels 
are created. We will explain only the latter here, since the first steps consist only of 
standard Mathematica transformations and are not special to CrasyDSE. Finally we 
explain some details for this specific example in the C++ code. The provided files 
allow the interested reader to follow the complete procedure, from the derivation of 
the DSEs to their numeric solution, in detail. 

5.1. Yang-Mills theory coupled to a scalar field 

In nature elementary matter fields are fermions. In quantum chromodynamics, 
for example, these are the quarks which interact via gluons. However, since their 
spin is 1/2, quarks are Dirac fields and consequently represented by spinors. An ad- 
vantage of functional methods is that they do not suffer from fundamental problems 
when dealing with anti-commuting fields. However, calculations are complicated 
by the Dirac structure, since it allows more dressing functions than for a simple 
scalar, see, e. g., [52] • While at the level of propagators this is still doable, see, 
e. g., [33ti36] . three-point functions become already quite tedious |37]. Since some 
non-perturbative phenomena like confinement may not depend on the fields being 
spinors or scalars, one can alleviate calculations by replacing the quarks by scalar 
fields. In order to mimic quarks such fields also have to be in the fundamental rep- 
resentation. The calculation used as an example for the presentation of CrasyDSE 
in this section is motivated by investigations along these lines [301 1381 - H2] . 

The renormalized action of this theory in Landau gauge reads in momentum 
space 





+ 



/ 



d qid q2 
(27r)2<i 



z Z, g r'^q,,Al{q{)A%q2)Al{-q, - q,) 



+ Z^gm2q2, + q^^)Al{q^)^\q2W{-q^- q^) +... , 



(10) 
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Figure 4: Top: The truncated DSE of the scalar two-point function. The full DSE can be found, 
for example, in ref. [38]. Bottom: The truncated DSE of the scalar-gauge field vertex. The second 
and third diagrams of the vertex DSE are called Abelian and non-Abelian diagrams, respectively. 
Gauge fields have red, continuous lines and scalar fields blue, dashed lines. Thick blobs denote 
dressed vertices. All internal lines are dressed propagators. 



where T^" are the Hermitian generators in the fundamental representation of SU {Nc) 
with the structure constants /"'"^. The dots correspond to four-point vertices and 
terms with Faddeev-Popov ghosts. The former are dropped in our truncation and 
the latter do not appear in the DSEs considered here, which are those of the scalar 
two-point function and of the scalar-gauge field vertex. The full DSE of the two- 
point function can be found, for example, in ref. [38]. Here we neglect for both DSEs 
all diagrams containing four-point functions which renders the scalar gap equation 
diagrammatically equal to the quark gap equation. This truncation is also motivated 
by the fact that two-loop diagrams are subleading in the UV. 

In the following we will focus on the scalar sector of the theory, viz. the scalar 
propagator and the scalar-gauge field vertex. As can be seen from the truncated 
DSEs of the scalar two-point function and the scalar-gauge field vertex in Fig. |4] 
we have the following four quantities left in our truncation: the propagators of the 
scalar and the gauge fields, the three-gauge field vertex and the scalar-gauge field 
vertex. 

The scalar propagator and the scalar-gauge field vertex are parametrized as 

iDs)^niP) = ^Smn, (11) 

where in the case of the vertex pi and p2 are the momenta of the scalar particles 
and z = pi ■ P2 / {\Pi\\P2\) ■ Introducing a sharp momentum cutoff A to regularize the 
self-energy contributions we need to approximate the three dressing functions 



A, : [0,A2] ^ M , (12) 
: [0, (A/2)']' X [-1,1] ^ M, (13) 

which will be done by linear interpolation for Ap^ p^ and linear interpolation as well 
as Chebyshev expansion in the case of Ag. Choosing the cutoff in the vertex to 
be smaller by a factor of two has the effect that only the dressing functions Ap^^p^ 
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will be called at large momenta outside their domain when evaluating the self- 
energy integrals. We will approximate them with their bare value Ap-^^p^ = Z\ where 
necessary. 

Additional information is required for the gauge field propagator and the three- 
gauge field vertex. For the latter we use for simplicity 

r^.r"^''=(Pi,P2,P3) = f r;^.f "^'•^'W(pi,p2,P3), (14) 

^3 

where finiteness of the ghost-gauge field vertex in Landau gauge Z\ = \ and the 
Slavnov- Taylor identity ZijZx = Z^/Z^ |13] have been used. The bare vertex 



^^uff''°''"^'^'^\pi, P2-, Ps) is given in Eq. (|20[). For the dressing function Z{p'^) of the 



Landau gauge field propagator we employ a fit to the solution of the ghost-gluon 
system obtained within the DSE framework provided in ref. [TT] (see also section 



-7 



Z{x) = {^] R'ix), (15) 



+ d 



QCD J \^^QCD 



2k 



2k ' 

l + c{^] +d' 



R{x) 

l + c(- 

where c = 1.269, d = 2.105, 7 = -13/22, k = 0.5953, Aqcd = 0.714 GeV and 



/ N tt(0) , , 



In 



with ai = 1.106, 03 = 2.324, h = 0.004, 62 = 3.169 and a(0) = 8.915/Nc, being 
the number of colors which we take to be three. Note that the renormalization 
constants Z3, Z3 can be calculated from the running coupling a as described in [TT]. 

For the present system the iteration procedure is very stable and we can start 
from a massless bare scalar propagator Ag = 1 and a bare scalar-gauge field vertex. 
As previously mentioned the integrals are regularized via an ultraviolet cutoff. To 
determine the renormalization constants of the scalar propagator we fix the values 
As(/i^) and Zmfn"^- Furthermore the vertex will be renormalized by enforcing the 
Slavnov- Taylor identity Z1/Z3 = Z^/Z^. With this prescription the system is then 
multiplicatively renormalizable, i. e., the expressions Z^Ag and {Zi)~^ Ap^^p^ are inde- 
pendent of the renormalization point /i^. Results confirming this for the propagator 
are shown in fig. [5j 

We can here also illustrate how extending truncations complicates the system of 
equations: A simple truncation takes into account only the scalar propagator. In this 
case we need the gauge field propagator and the scalar-gauge field vertex as input 
and we only have to calculate one integral. Its integrand is comparatively simple. 
Going only one step further by including also the scalar-gauge field vertex requires 
solving in total for three dressing functions (the vertex has two and the propagator 
one) by calculating five integrals, whereby the complexity of the integrands has also 
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Figure 5: Dressing function of the propagator times its renormalization constant. Multiplicative 
renormalizability of the system propagator and vertex is clearly visible as the results obtained with 
two different choices of /i^ are indistinguishable. Independence from the approximation method is 
seen from the results obtained with the Chebyshev expansion method. For comparison the result 
with a bare scalar-gauge field vertex is included. 



increased considerably. 

We should mention that the employed truncation is not state of the art but 
it is sufficient to illustrate many features of CrasyDSE. More elaborate results for 
Yang-Mills theory coupled to a fundamental scalar will be presented elsewhere [30] . 

5.2. Generating the C++ files with the integrands 

We turn now to the creation of the kernels with Mathematica. The subsequent 
explanations follow the notebook CrasyDSE.YM+scalar.nb. To initialize the re- 
quired functions and the expressions for the DSEs we evaluate the initialization 
cells with 

FrontEndTokenExecute ["Evaluatelnitialization"] 

or via the menu entry Evaluation -+ Evaluate Initialization cells. Now the vari- 
ables containing the expressions of the integrals are defined: gapAlgProj Loop- 
Integrand is the integral of the gap equation and vertexAbelianProjplFinal, 
vertexNonAbelianProjplFinal, vertexAbelianProjp2Final and vertexNonAbe- 
lianProjp2Final correspond to the integrals of the Abelian and non-Abelian di- 
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agrams projected onto the external momenta pi and p2- Furthermore the package 
CrasyDSE is loaded. 

Basically for the generation of the C++ files only one function is needed. How- 
ever, before wc can use it we need to define several expressions. They are split into 
lists containing parameters, momentum variables and dressing function names. In 
the following the order of the elements in lists is very important. Thus one has to 
be very careful when one changes something later on, because this could lead to 
inconsistencies with the C++ code. 

In the gap equation two parameters appear: the number of colors Nf. and the 
coupling constant g. We put them both into a list: 

parasGap = {Nc, g}; 

Furthermore we specify the external (ps:= p^) and internal (qs:= (f, ct:— cos((/7) = 
p ■ q/\p\\q\) variable names: 

extVarsGap = {ps}; 
intVarsGap = {qs, ct}; 

Although we did not introduce abbreviations of momentum combinations for the 
gap equation, we need to define a variable for this, because we will need it later: 

extraVarsCListGap = {}; 

Finally we have to specify which dressings appear in each equation. They are sep- 
arated into two lists, one for the dressings of the Green function we are calculating 
and one for all other dressings belonging to other Green functions. The propagator 
has only one dressing function Dg, so the first list contains one item only: 

dressingsGap = {Ds}; 

But the gap equation also depends on other dressing functions, namely on the one of 
the gauge field, Da^ and on the two of the scalar-gauge field vertex, d'^^I- and D^^l-. 
All dressings belonging to the same Green function have to be grouped together in 
one sublist: 

otherGreenFuncsGap = {{DA}, {DAssbl, DAssb2}}; 

These are the lists required for the gap equation as input for generating the C++ 
files. 

The lists for the vertex equation have the same structure and we only list them 
here: 

parasVertex = {Nc, g}; 
extVarsVertex = {pis, p2s, ca}; 
intVarsVertex = {qs, ctl, ct2}; 

extraVarsCListVertex = {{plp2, ca Sqrt[pls p2s]}, 
{plq, Sqrt[pls qs] Cos[ct2]}, 

{p2q, (ca ct2 + Sqrt[l - ca~2] ctl Sqrt[l - ct2"2] ) Sqrt [p2s qs]}, 

{plmqs, pis - 2 plq + qs}, 

{p2mqs, p2s - 2 p2q + qs}, 

{plmp2s, pis + p2s - 2 Sqrt [pis p2s] ca}}; 
dressingsVertex = {DAssbl, DAssb2}; 
otherGreenFuncsVertex = {{DA}, {Ds}, {DAAA}}; 
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Note that the vertex has two dressings by itself and depends on three other Green 
functions. Furthermore wc provided with the hst extraVarsCListVertex the defi- 
nitions of employed abbreviations, e. g., plp2:= pi ■ P2 = cos{a)\pi\\p2\. 

Before we generate the C++ files we split off numeric coefficients of the inte- 
grands. We call the resulting expressions kernels and coefficients. Instead of doing 
this by hand, we use the function split Integrand. It takes as arguments an expres- 
sion and a list of variables. Everything in the overall factor that does not contain a 
variable will be put into the coefficient: 

{coeffGap, kernelGap} = 
split Integrand [gapAIgPro j Looplntegrand , 

Join[extVarsGap, extraVarsCListGap[[All, 1]], intVarsGap]] /. 
Zlh :> 1 // Simplify 

— > {(g~2 (-1 + Nc~2))/( 

8 Nc \[Pi]"3), (l/pplusqs)(l - ct-2)-(3/2) 
DA[qs] (DAssblEps, 

pplusqs, -((ps + ct Sqrt [ps qs] ) /Sqrt [pplusqs ps])] - 
DAssb2 [ps, 

pplusqs, -((ps + ct Sqrt [ps qs] ) /Sqrt [pplusqs ps])]) Ds [pplusqs]} 

We discarded the renormalization function Zlh here since its implementation is 
handled manually. 

For the vertex we do the same but bear in mind the following structure: Both for 
coefficients and kernels every loop integral is treated as a single expression, and for 
every equation all loops are grouped into lists. The syntax of the list of coefficients 
or kernels is thus 

{{loop 1 of eq. 1, loop 2 of eq. 1, ...}, 
{loop 1 of eq. 2, loop 2 of eq. 2, ...}, ...} 

For the first and second projections we split the integrands as follows: 

{coef f sVertexProjpl , kernelsVertexProjpl} = Transpose [ 
splitlntegrand [# , 
Join[extVarsVertex, extraVarsCListVertex [ [All , 1]], 
intVarsVertex] ] & /@ {vertexAbelianProjplFinal , 
vertexNonAbelianProjplFinal} /. {Zlh :> 1, Zl :> 1} // Simplify]; 
{coef f sVertexProjp2, kernelsVertexProjp2} = 
Transpose [ 
splitlntegrand [#, 
Join [extVarsVertex, extraVarsCListVertex [ [All , 1]], 
intVarsVertex]] & /@ {vertexAbelianProjp2Final, 
vertexNonAbelianProjp2Final} /. {Zlh :> 1, Zl :> 1} // Simplify]; 

Again we have discarded the renormalization functions. The final lists of kernels 
and coefficients are 

kernelsVertex = {kernelsVertexProjpl, kernelsVertexProjp2}; 
coeffsVertex = {coef f sVertexProjpl , coef f sVertexProjp2} 

--> {{g/(16 Nc Pi-3), -((g Nc)/(32 Pi-3))}, 
{g/(16 Nc Pi-3), -((g Nc)/(32 Pi-3))}} 
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We show the coefficients exphcitly. One can easily spot the and Nc dependences 
of the Abehan and non-Abehan diagrams, respectively. 

Finally we have everything to generate the C+ + code. We do so with the function 
exportKernels: 

exportKernels [{FileNameJoin[{NotebookDirectory [] , ".."}], 
"kernelsAll"}, 
{"scalar.QCD.hpp"}, 
{{{"coeffGap", "kernelGap"}, 

{coef fGap}, 

{kernelGap} , 

dressingsGap, 

otherGreenFuncsGap , 

parasGap , 

extVarsGap , 

intVarsGap , 

extraVarsCListGap} , 
{{"coef f sVertex" , "kernelsVertex"} , 

coef f sVertex, 

kernelsVertex, 

dressings Vertex , 

otherGreenFuncsVertex , 

parasVertex, 

extVarsVertex , 

intVarsVertex , 

extraVarsCListVertex} 
}] 

It will create two files kernels All. hpp and kernels All. cpp. The filenames are deter- 
mined by the first argument where we also indicate that the files should be exported 
to the parent directory. The second argument here is the name of an additional 
header file which contains functions specific to this example. The third argument 
contains all the information we gathered above: It is a list where every item corre- 
sponds to one DSE. For every DSE we have the following entries: 

• The C++ names of the functions containing the coefficients and the kernels. 

• A fist with the expression (s) for the coefficient (s). 

• A fist with the expression(s) for the kernel(s). 

• The list of dressings for this Green function. 

• The list of dressings from other Green functions. 

• The list of parameters. 

• The list of external variables. 

• The list of internal variables. 

• The list of extra variables. 
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Note that without specifying a path in the first argument of exportKernels the 
files will be created in the directory of the notebook. 

We want to mention here the function f unctionToString, which is used by 
exportKernels but can also be used directly by the user. It creates a string of the 
expression given as its argument similar to the Mathematica function CForm, but it 
replaces some common functions like Power or Sin by its C++ counterparts pow or 
sin. If a function is not included, it can be added by hand, for example: 

f unctionToString [a~b + Sin[b]/10 - 5 Sinh[a] , {Sinh :> sinh}] 

— > 0.1*(sin(b)) + -5.*(sinh(a)) + (pow(a, b)) 

This finishes our work in Mathematica and we proceed with the C++ code. 

5.3. Numerical code 

The C++ code, contained in scalar_main.cpp, is extensively commented and 
every variable that appears is described directly in the file. Here we only give a 
rough overview of the required initializations. 

In the file scalar^main.cpp first the model and then all Green functions are ini- 
tialized. For the former the definitions are as simple as providing numeric values for 
some parameters, e. g., A^^^ = 3. All Green functions are defined as a dse structure, 
which contains a pointer to mod which is reserved for hosting model parameters, 
see also Fig. [3j Since the gauge field propagator and the three-gauge field vertex 
are given by ansatze the only additional information these structures need are the 
corresponding definitions. The dynamically calculated Green functions also contain 
an array of quad structures, namely one for each different integration. Consequently 
variables like the numbers of integration points and the quadrature types have to 
be initialized. We also have to provide starting expressions for the dressings and in- 
formation on the other Green functions contained in a DSE - handled via the array 
otherGF. For example, for the gap equation these are the gauge field propagator 
and the scalar-gauge field vertex. In the C++ code dressing functions do not have 
a specific name, but are just collected in the function dress where it is important 
at all times to maintain the same assignment of the dressing functions as in the 
notebook. The integrands of the self-energy are defined in the kernels file created 
with CrasyDSE^YM+scalar.nb. Also a renormalization procedure has to be defined. 
Finally, the iteration is done with the function meta_solve_iter. 

6. Landau gauge Yang-Mills theory 

As a second example we use pure Yang-Mills theory which requires some different 
methods. The most obvious change is that we use a Newton's method instead of a 
direct fixed point iteration. Again we provide the complete Mathematica and C++ 
code together with the program. However, as the creation of the kernel files is rather 
similar to the case of the previous section we refrain from showing any details. 

6.1. Truncation and ansatze 

As in the last section we will employ the Landau gauge. The DSE system 
truncated at the level of propagators has been investigated with DSEs for some 
time now, see, for example, [TUl [T21 El EZl EH Hi]- We will here reproduce the 
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Figure 6: The truncated two-point DSEs of pure Yang-Mills theory. Gluons have red, continuous 
lines and ghost fields green, dashed lines. Thick blobs denote dressed vertices. All internal lines 
are dressed propagators. 



solutions of refs. EH HS]. Besides employing a Newton procedure to solve this 
set of equations another difference to the previous section lies in the renormalization 
procedure: Here we work with subtracted DSEs. 

The system we investigate consists of the ghost and gluon two-point DSEs. The 
former is used without change, while the gluon DSE is truncated [101111]: We neglect 
all diagrams involving a bare four-gluon vertex, i.e., the tadpole diagram and all two- 
loop diagrams. They are subleading in the UV and it was shown analytically for 
the scaling solution that they are also subleading in the IR 06]. The remaining 
unknown quantities are the ghost-gluon and the three-gluon vertices for which we 
use suitable ansatze. The truncated set of DSEs is depicted in Fig. |6j 

The ghost and gluon propagators are given by 

Kip) = {a,. - ^) (17) 

G^\p) = -S-b9ifl, (18) 

For the ghost-gluon vertex T^'^'^''^'"'{pi, p2, ps) the bare version, 

rf-'-^'^''~'\p,,p,,p,)=zgr'^p,„ (19) 

is used as motivated originally by an argument of Taylor. However, several stud- 
ies both on the lattice 07] and in the continuum [IHl SH] confirmed this to be 
a very reliable ansatz. For the full three-gluon vertex r^^^'"^'^(pi,p2,P3) we use 

the tensor structure of the bare vertex T'i]iutf''"''"^'^^\pi,P2,P3) amended by a dress- 
ing D^^^{p1,P2,P2) that guarantees the correct UV behavior of the gluon dressing 
function [31] : 

j.AAA,abc,m ^p^^p^^ p^^^^g jabc ^^^^(^^ - pi) ^ + g,^{p, - p2) , + g,,{p, - pa).) , 

(20) 

r^.r"''"^(pi'P2,P3) = r;:i''''^''^'(°)(pi,P2,P3)/^^^^(j>Lj>LpD, (21) 

^1 [Z[pi)Z[pi)) 
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Figure 7: Schematic illustration of structures defined in the code and their dependencies for the 
pure Yang-Mills system. 



a is a parameter chosen as 36, where 6 is the anomalous dimension of the ghost 
propagator. Zi is the renormahzation constant of the three-gluon vertex. The 
dependencies of all Green functions on each other are shown in Fig. [7} 

An important issue of the gluon DSE are spurious quadratic divergences. They 
appear because we employ a numerical cutoff as UV regularization which breaks 
gauge invariance. There are several ways to deal with them, see, for example, 
[T^ 1211 ED]- Here we subtract an additional term in the kernel of the gluon loop in 
the gluon DSE [31]. The derivation of the DSEs with DoFun and the projection to 
scalar quantities are described in the notebook DoFun^YM^^d-nh and the creation 
of the kernel files in CrasyDSE^YM^^d-nh. For details we refer to them. 

6.2. Renormalization and solution 

For the present system we will use subtracted DSEs, i. e., we subtract from a 
DSE at external momentum p the DSE at a fixed external momentum p^: 

D-'{/) = Z-' + Tl{p') D-\p') = D'\pI) + TI{p')-TI{pI). (23) 

Thus we can trade the renormalization constant Z for specifying the value of the 
dressing function at the subtraction point po- For the gluon propagator we choose 
the subtraction point po at sufficiently high momenta since we expect the two-point 
function to be divergent at low momenta. For the ghost, however, it is most ad- 
vantageous to specify the dressing at zero momentum. These conditions are bound- 
ary conditions for the integral equations. As it turns out two different types of 
solutions can be found depending on the value of the ghost dressing at zero mo- 
mentum: Choosing finite values a solution of the family of decoupling solutions 
emerges HU EI] , while with an infinite zero-momentum dressing we get the scal- 
ing solution [ini 12]. The former has a finite gluon propagator and a finite ghost 
dressing function at zero momentum and the latter an IR vanishing gluon propa- 
gator and an IR divergent ghost dressing function. Thereby the divergence of the 
ghost dressing and the vanishing of the gluon dressing can be described by power 
laws whose exponents 5gh and 5gi, respectively, are related by 5gi + 25gh = 0, whereby 
Sgh := K = 0.595353 can be calculated analytically [IHl 152] . 
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There are several choices at which point of the calculation the subtraction of a 
DSE can be performed. For illustration purposes we employ a different one for each 
propagator: For the ghost we use the subtracted expression in the kernel file. This is 
advantageous because the limit of vanishing external momentum can be done analyt- 
ically but is problematic numerically. For the gluon propagator, on the other hand, 
we only create the unsubtracted expressions. The subtraction is then performed 
with a function in C++. Here the subtraction is numerically unproblematic. 

The specific renormalization procedure has to be taken into account in the renor- 
malization function in the C++ code. Furthermore, it is important to note that this 
function does not calculate the right-hand side of a DSE as in the example of the 
scalar system but the difference between the right- and left-hand side of a (sub- 
tracted) DSE: 

E{p') := -D-\p') + D-\pI) + n(p2) - n{pl). (24) 



The reason is the employed Newton procedure as described in section |4.3.2| which 
attempts to bring -E(p^) to zero. E{p^) is calculated for every external momentum 
and saved as an array of the DSE structure. 

The behavior of the dressing functions in the IR and UV is known analytically. In 
the intermediate regime, between two given momenta e and A, above and below the 
IR and UV cutoffs, respectively, they are expressed by an expansion in N Chebyshev 
polynomial^ 



N-l 



Gjm{p') = exp J2 ^ UM{p')), (25) 

1=0 
N-l 

Zjm{p') = exp J2 T,{M{p')), (26) 



1=0 



where M{p'^) maps the regime [e. A] to [—1, 1]. For momenta below e we employ a 
power law with the given exponent and the coefficient calculated from the lowest 
known point in the Chebyshev expansion: 

GiR{p^) = A^^^\py^\ (27) 
Z,^(/) = AW(p2)5,,_ ^28) 

For momenta higher than A an extrapolation in agreement with the UV behavior is 
chosen: 

Guvip") = G{s'){w\og{p'/s') + 1)^ (29) 
Zuvip") = Zis'){w\ogipys') + ly. (30) 

s is the highest momentum at which the Chebyshev expansion is known, 6 or 
7 are the anomalous dimensions of the ghost and gluon, respectively, and w = 
11 Nca{s) G{sY Z{s)/12n . a(p^) is a possible non-perturbative definition of the 



®The exponential is chosen due to better convergence properties. For such an expansion see 
also, for example, ref. [551 [55]. 
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Figure 8: Ghost (left) and gluon (right) dressing functions. The continuous hne corresponds to 
the scahng solution, the dashed lines are solutions of the decoupling type. 



running coupling 



(31) 



The value of Q;(/i) is an input parameter and sets the scale: At /i we have G'(/i^)^Z(yU^) 
1. 

As starting functions for the propagator dressings in the intermediate regime we 

use 



Zansip') = flRipTipT" + CUV " fuv{p'), 
Gansip')=flRip')ipy^' + l, 

with the IR and UV damping factors given by 

LiR 



fmip') 
fuvip'] 



LiR + p'- 

.2 



Luv + p"^ 



(32) 
(33) 



(34) 
(35) 



where L/^j and Ljjy are dimensionful parameters conveniently set to 1. The parame- 
ter cuv can be used to adjust the starting function in order to speed up convergence. 
Here it is chosen as 1. These ansatze respect the qualitative IR behavior which leads 
to a faster convergence than starting with constant functions. In general the Newton 
procedure becomes more stable when using starting functions close to the solution. 
The starting functions being not differentiable at the UV matching poses no prob- 
lem. 

For the integration we used Gauss-Legendre quadratures. Furthermore it was 
advantageous to split the radial integration at the value of the external momentum 
which requires setting bound_type of all quadratures to 1. Note that this has to be 
done after init_quad, which sets the default value 0. This allows a higher precision 
of the final result. However, this comes at the prize of slowing down the integration 
as the integration boundaries have to be calculated for every external momentum. 

In Fig. [8] we show the results of the calculations for different boundary conditions 
of the ghost. It is clearly visible that all solutions coincide in the UV and only show 
their distinct behavior in the IR. In table |3] we provide the input parameters used 
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parameter 


1 

value 


parameter 


1 

value 


Nc 


3 




0.595353 




1, 0.5 


h 


10-^ 


UV cutoff 


10^ 


e 


2 X 10-^ 


IR cutoff 


10-12 


A 


0.99 X 10^ 


6 


-9/44 


gluon subtraction point po 


1.2 


7 


-13/44 


value of gluon dressing at po 


0.93 


LiR 


1 


value of ghost dressing at p = 


0, 5, 10, 25 


Luv 


1 







Table 3: Parameters for the calculation. Where several values are given see text for details. 




Figure 9: Left: Ghost and gluon propagators. The continuous (red) lines corresponds to the 
gluon, the dashed (green) lines to the ghost. Thick and thin lines corresponds to dilferent choices 
of a{ix^). Right: Coupling for the two choices of a(^^). The two lines (black straight and gray 
dashed) are almost indistinguishable. 

for our calculations. The reached precision can be seen from how well the results 
fulfill the DSEs, i. e., how close Eip^) approaches zero. Its norm goes with double 
precision down to about 10-^. Using long double variables instead this value can 
be made even lower. 

We also tried a second choice for a{jjp') to test the code, namely a{ji^) = 0.5. As 
expected the propagators change by a constant factor due to multiplicative renormal- 
izability, whereas the running coupling a{p^) is independent of /i^, see, for example, 
[T2l [55] . A comparison between a{n'^) = 1 and 0.5 is depicted in Fig. [9] 

Finally we want to make some technical remarks: With these examples we only 
want to illustrate the basic use of CrasyDSE so the code is not optimized and we 
expect that the runtime can be improved considerably. Another point is that we 
also tried a simple fixed point iteration but did not get a solution. This may indicate 
that this method is not suited for this problem. 

7. Summary and outlook 

The strength of the framework provided by CrasyDSE lies in its ability to handle 
large numbers of dressing functions. Thus one of its fields of application is the exten- 
sion of current truncation schemes by including higher vertices and/or enlarging the 
tensor bases of Green functions. For example, it is expected that going beyond the 
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lliljC-i jJvjiciLiwiio . 


1 iTi fifi y 

llliCdl 


expansions: 


Chebyshev 


solution methods: 


fixed point iteration, Newton 


quadratures: 


Gauss-Legendre, Chebyshev, Fejer, double exponential 



Table 4: Currently implemented numerical methods. 



current truncation schemes in the Landau gauge makes DSE solutions more compet- 
itive to lattice solutions in the mid-momentum regime. Since the number of dressing 
functions increases at non-zero temperature and/or non-zero density corresponding 
calculations can profit from CrasyDSE too. Finally there are also interesting cases 
for which no numerical calculations have been done successfully yet because their 
systems of DSEs are very complex. 

Combining CrasyDSE with DoFun there exists a sound framework for the treat- 
ment of all such systems of DSEs, from their derivations to their numeric solutions. 
Furthermore, we plan to extend CrasyDSE by adding further approximation meth- 
ods or new solving algorithms to the ones given in table |4] as required by individual 
cases. 
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Appendix A. Initializing and accessing x A and A 

Here we provide some details on accessing and initializing the arrays double 
*x_A and double *A which are members of the structure struct dse. Details on 
their structure can be found as comments in DSE.cpp. Here we only describe the 
functions to access them properly which should be sufficient for most applications. 
For every DSE we also have to define the number of external momenta and the 
number of dressing functions, int dim_x and int dim_A, respectively. For every 
external momentum the number of interpolation points or expansion coefficients 
has to be specified. The corresponding values form the array int *n_A. As an 
introductory example for the use of these variables one can consider the interpolation 
example interp^main. cpp. 

Appendix A.l. The array x_A 

The array x_A contains the interpolation points and/or Chebyshev nodes de- 
pending on the chosen interpolation method. In the case of Chebyshev interpola- 
tion only the discrete variables need to be initialized by hand in x_A. However when 
using the DSE solving routines it is necessary that also the continuous Chebyshev 



26 



nodes are stored in x_A. These are automatically initialized in x_A by calling the 
function void Cheb_init_cont_xA(void *dse_param). In the provided examples 
Cheb_init_cont_xA is called in the member void (*init_dress) of the structure 
dse , e. g., scalprop_init_dress_cheb for the scalar propagator. Note that due 
to technical reasons the values for each continuous variable are ordered from low to 
high for a linear interpolation, but from high to low for a Chebyshev expansion. 

Independent of the interpolation method the discrete external momenta have to 
be initialized by the user in x_A. In the case of linear interpolation this is also true for 
the continuous argument^ Assume we are interested in the grid point i. Its grid 
coordinates are stored in the array *ind by calling the function void index (int 
*ind, int *n, int dim, int i ), where *n are the number of grid points in every 
external momentum, given by n_A, and dim is the number of external momenta, given 
by dim_x. Then we can obtain the index of the idir-th component of the i-th grid 
point with the function int xA_index(int *ind, int idir, void *dse_param) . 
When x_A is correctly initialized access to the coordinates of the i-th grid point is 
provided via the function void outer .argument (int i, void *dse_param) which 
stores them in outer.arguments. 

Appendix A. 2. The array A for linear interpolation 

When linear interpolation is used the array A contains the function values at the 
external momenta x_A. The function int A_index(int *ind, int idress, void 
*dse_param) returns the index of the grid point with the coordinates ind obtained 



with the help of the function index, see section [Appendix A.l[ The argument 
idress is the index of the dressing function. 

Appendix A. 3. The array A for Chebyshev expansion 

For a Chebyshev expansion the array A contains the Chebyshev coefficients. They 
can be initialized from a (user-provided) function f unc by calling void Cheb_init- 
_coeff_mult (double (*f unc) (double *x, void *param) , void (*traf) (double 
*x, void *paraiii) , int dim, int *n_f, int *n_c, double *c_ar, void*param). 
The argument traf defines the transformation of the domain of the Chebyshev 
polynomials ([—1, 1] and direct products thereof) to whatever is the domain of 
the function. We recommend the function void Cheb_gen_traf (double *x, void 
*param) for this purpose which transforms the interval depending on int *cheb- 
_trafo, a member of the corresponding dse structure. If cheb_traf o [i] is 0/1/2 
no transformation/a linear transformation/a logarithmic transformation is employed 
for the i-th external momentum. In what follows we give the expressions usually 
used as arguments of Cheb_init_coef f jnult in parenthesis, dim defines the num- 
ber of continuous arguments (dim_x - dimjnat), *n_f (n_A) is the number of points 
used for evaluating the inner product to calculate the *n_c (n_A) Chebyshev coeffi- 
cients stored in *c_ar (A) and param are some parameters (dse_param). In general 
n_f=n_c should be used. We note here that if more than one dressing function 
is to be interpolated the Chebyshev coefficients can be stored in one array A by 
simply passing the argument A+idress*ntot_A as *c_ar where idress denotes the 
idress-th function. 



^For Chebyshev interpolation this has to be done manuaUy as well if one wants other external 
momenta than the Chebyshev nodes initialized by Cheb_init_cont_xA. 
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Appendix B. Structure of user-defined functions 

For the following members of dse or quad structures the user has to provide 
functions: 

• void def .domain (double *x, void *dse_paraiii) : Here the user defines the 
boundaries of the interpolation domain for the dressing functions. The bounds 
are saved in the array double *domain, a member of the structure dse_param. 
The lower and upper bounds of the i-th external momentum are saved in 
domain [2*i] and domain [2*i+l] , respectively. The array x refers to the ex- 
ternal variables. As examples consider the following functions: scalgluevert- 
_def .domain in scalar.QCD.cpp or interp.def .domain in interp.cpp. def- 
_domain is a member of a dse structure. 

• double interp_of f domain(int *pos, double *x, int iA, void *dse_param) : 

This function is required for extrapolation, i. e., when the domain of the 
interpolation of a dressing function is left. The array pos of length dim_x 
contains the values 0, 1 and 2 for every external momentum. If pos [i] =0, 
the interpolation domain for the i-th external momentum is not left, while 
a value of 1 or 2 means that the lower or upper bounds are crossed, re- 
spectively, pos is created automatically based on the information provided 
in def_domain. iA tells which of the dim_A dressing functions is required 
and the array x contains the values of the external momenta. dse_param 
refers to the dse structure to which this dressing function belongs. Ex- 
amples are the functions ghprop_interp_of f domain_cheb in YM^d.cpp or 
scalvert_interp_off domain in scalar.QCD. cpp. interp_off domain is a mem- 
ber of a dse structure. 

• void (*renorm) (void *dse_param): This function has to perform several 
tasks. It does not only implement the renormalization procedure but also 
contains all steps required to proceed from the results of the integration to the 
expressions required by the solving algorithms. 

The self-energy contributions of single diagrams are calculated by the solv- 
ing algorithms with the function self energy. It stores the results for each 
dressing function and diagram in the dse member double *self_A. This 
array is organized as follows: For every quadrature Q [i] a block of size 
n_loop [i] *dim_A*ntot_A is used. These blocks are separated into ntot_A 
sub-blocks of size dim_A*n_loop [i] , i. e., each entry contains the result of the 
integration of one diagram for a specific external momentum. Appropriately 
summed up we obtain the result for the right-hand side of a DSE for all ex- 
ternal momenta. This summation as well as the renormalization have to be 
done in renorm. Furthermore, if the DSE is not projected directly onto its 
dressing functions, the linear system of equations to obtain them has to be 
solved here. Finally, depending on the solving algorithm for the DSEs, either 
the dse members A or E have to be set: In the case of a fixed point iteration the 
results can be directly saved to A if a linear interpolation is used. For a Cheby- 
shev expansion the new coefiicients A are calculated with Cheb_coef f _mult. If 
the solving algorithm is Newton's method the renormalization function calcu- 
lates of Eq. (pj). Examples are scalgluevert_renorm in scalar_QCD.cpp 
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and prop_renorni_cheb_secant in YM^d.cpp. renorm is a member of a dse 
structure. 



void integrand (double *erg, double *x, void *int_param): The actual 
kernels for the integration are contained in this function. Every quad struc- 
ture handles the integration of one or several integrals and the values of the 
integrands at the integration variables given by the array x are stored in the 
array erg. The external momenta are accessed via int_param, which refers 
to a dse structure. Its member double *outer_arguments has to contain 
the external momenta, integrand can be a user-written function, but more 
commonly it will be created by the Mathematica functions of CrasyDSE. Ex- 
amples are sphere_integrand in sphere. cpp, which was created manually and 
is thus relatively simple, and kernelsVertex in kernels All. cpp of the scalar 
QCD example, which was created in CrasyDSE_YM+scalar.nb. integrand is 
a member of a quad structure. 

double Jacob (double *x) : This function can be used for the Jacobian of the 
integral measure. Since the Jacobian is automatically taken into account for 
every integration, it must always be defined. Thus, if it is already contained in 
the integrand, this function must return 1. The array x holds the integration 
variables. Examples are sphere_jacob in sphere. cpp and ghprop_jacob in 
YM^d.cpp. Jacob is a member of a quad structure. 

void coeff (double *coef f icients , void *int_param): Any trivial coef- 
ficients of the integrals which do not depend on any momenta can be put 
into this function. Their values are saved in the array coefficients which is 
multiplied with the results from the integration. int_param refers to the quad 
structure of the integration. Similar to integrand this function can be writ- 
ten manually or created with Mathematica. Examples are sphere_coeff in 
sphere. cpp, which was created manually, and coeff sVertex in kernels All. cpp 
of the scalar QCD example, which was created in CrasyDSE_YM+scalar.nb. 
coeff is a member of a quad structure. 

void boundary (double *bound, double *x, int idim, void *int_param): 
The limits of the integration are defined in this function. For the idim-th 
integration variable bound [0] and bound [1] are set to the lower and upper 
integration bounds, respectively. The array x contains the external momenta, 
on which the integration bounds may depend. If this is the case, the variable 



bound_type has to be set to 1, see section 4.2 int_param refers to the quad 
structure of the integration. Examples are sphere_boundary in sphere, cpp 
and scalgluevert .boundary in scalar. QCD. cpp. boundary is a member of a 
quad structure. 

Appendix C. Overview of all C++ functions 

This appendix provides tables with all C++ functions of CrasyDSE relevant for 



the user. Table C.5 contains the most prominent functions of CrasyDSE, table C.6 



the functions and variables for approximating the dressing functions, table C.7 the 



functions and variables for the integration and table C.8 further functions and vari- 



ables of dse structures. Details on the functions' arguments are provided in the 
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code in the function descriptions. For the user's convenience we provide the file 
template -DSE.cpp in the directory main which contains a hst of all functions and 
variables that have to be defined for a dse structure. 

Tabic C.5: Overview of the main C++ functions of 
CrasyDSE. 



Function 


File 


Short description 


int 


DSE.cpp 


index of coefficient with co- 


A index 




ordinatrs ind of dress-th 

\j 1. V_XXXX<_Ai U V-ykJ _L. XX Kj L \. m t \^ UX± 


(int *ind, int dress, void 




dressing function 


*dse_parain) 






void 


dressing, cpp 


Chebyshev coefficients from 


Cheb_coef f jnult 




array of function values at 


(int dim, int *n_f, int 




Chebyshev nodes 


*n_c , double *f_a.r, double 






*c_ar) 






double 


dressing, cpp 


interpolating Chebyshev 


Cheb gen dress int erp 




polynomial 


(double *x, int i, void 






*dse Daram^ 






void 


dressing, cpp 


maps domain of Chebyshev 


Cheb_gen_traf 




polynomials to domain of 


(double *x void *Daraiii^ 




function 


void 


dressing, cpp 


Chebyshev coefficients from 


Cheb_init_coef f jiult 




function 


(double (*funcl (double 






*x, void *paraiii) , void 






(*traf) (double *x, void 






♦param) , int dim, int *n_f. 






int *n_c, double *c_ar, void 






*param) 






void 


dressing, cvv 


initializes the continuous 


Cheb_init_cont_xA 




external variables to the 


(void *dse param) 




Chebyshev nodes 


void 


DSE.cpp 


deallocate members of dse 


dealloc_A_xA 




structure 


(void *dse_param) 






void 


quadrature, cpp 


deallocate members of quad 


dealloc.quad 




structure 


(void *quad_param) 






void 


quadrature, cpp 


transformed nodes and 


get_nw 




weights 


(double *x, double *w, 






int i, int idim, void 






♦bound params , void 






*quad_param) 






Function 


File 


Short description 



30 



Function 


File 


Short description 


void 


function, cpp 


coordinates of index- 


index 




th entry of a one- 


(int *i, int *n, int dim, 




dimensional array in a 


int index) 




n[0] X . . . xn[dim-l] grid 


void 


DSE.cpp 


allocate members of dse 


init_A_xA 




structure 


(void *dse_parain) 






void 


DSE. cpp 


default values of parameters 


init_dse_def ault 




for allocation of dse struc- 


(void *dse_parain) 




ture 


void 


quadrature, cpp 


allocate members of quad 


init_quad 




structure 


(void *quad_param) 






void 


quadrature, cpp 


integrate a set of functions 


integrate 






(double *erg, void 






*quad_parain, void 






*int_parain) 






double 


dressing, cpp 


linear interpolation of array 


L in_gen_dre s s _int erp 




of function values 


(double *x, int i, void 






*dse_parain) 






void 


DSE.cpp 


solves coupled set of DSEs 


meta_solve_iter 




by iteration 


(struct dse *DSE, int ndse. 






double epsabsstop, double 






epsrelstop, int iterstop, 






int output) 






void 


DSE. cpp 


sets outer_arguments to 


outer_argument 




the i-th external momenta 


(int i, void *dse param) 






void 


DSE.cpp 


prints version and author 


plot_credits() 




information 


void 


dressing, cpp 


output of dressing functions 


plot_dressing 




and interpolation points on 


(void *dse param) 




screen 


void 


DSE.cpp 


solves one DSE by iteration 


solve_iter 






(void *dse_param, int 






output) 






void 


DSE. cpp 


solves array of DSEs with 


solve_iter_secant 




Newton's method 


(struct dse *DSE, int ndse, 






int maxiter, double eps_E, 






int output) 






Function 


File 


Short description 
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Function 


File 


Short description 


void 


dressing, cpp 


output of dressing functions 


write_dressing 




and interpolation points in 


(void *dse_paraiii, char 




file 


*name, double *addpara, 






int length) 






int 


DSE. cpp 


index of idir-th direction 


xA_index 




of a grid point with coordi- 


(int *ind, int idir, void 




nates ind 


*dse_param) 






Function 


File 


Short description 



Table C.6: User defined members of a dse structure rel- 
evant for the approximation of functions. References to 
the scalar-gauge field vertex refer to the example of sec- 
tion [5} those to the ghost-gluon system to the example 
of section [H 



Member 


General purpose 


DSE usage/example 


double 
*A 


interpolation coefficients 


e.g., coefficients of Ap^ and 
for the scalar-gauge 
field vertex 


int 

*cheb_f unc_traf 


defines transformation of 
approximated function for 
Chebyshev expansion; 0: 
none, 1: logarithmic 


e.g., approximation of the 
exponential of the ghost 
dressing function instead of 
the dressing function itself 


int 

*cheb_traf 


defines transformation func- 
tion to interval [—1, 1] for 
Chebyshev expansion; 0: 
none, 1: linear, 2: logarith- 
mic 


e.g., logarithmic transfor- 
mation of definition interval 
of ghost dressing function 


double 
*cheb_x 


necessary for Chebyshev in- 
terpolation with diimiat> 



see interp.cpp for correct us- 
age 


void 

(*def _domain) 
(double *x, void 
*dse_paraiii) 


defines domain of definition 


e.g., [0,(A/2)^]' X [-1,1] 
for scalar-gauge field vertex 


int 
dini_A 


number of interpolated 
functions 


number of dressing func- 
tions for a given Green func- 
tion, e.g., 2 (Apj and Ap^) 
for the scalar-gauge field 
vertex 


int 

dimjiat 


number of discrete argu- 
ments of interpolated func- 
tions 


number of independent 
Matsubara frequencies 


Member 


General purpose 


DSE usage/example 
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Member 



General purpose 



DSE usage/example 



number of arguments of 
dressing functions, e.g., 3 
{pI,P2 and z) for the scalar- 
gauge field vertex 



int 
dini_x 



number of arguments of in- 
terpolated functions 



double 
(*dress) 

(double *x, int i, 
void *dse_paraiii) 



returns value of the i-th 
dressing function at the mo- 
menta X 



e.g., Ap^ and Ap^ for the 
scalar-gauge field vertex 



double 
*E 



array of the values on the 
left-hand side of Eq. ^ for 
Newton's method 



has to be calculated in the 
renormalization procedure 
defined by the user, e.g., in 
prop_renorm_cheb_secant 
for the ghost-gluon system 



void 

(*init_f unc) 
(double *x, void 
*dse_paraiii) 



double 

(*interp_off domain) 
(int *pos, double 
*x, int iA, void 
*dse_param) 



function used for initializa- 
tion of the Chebyshev coef- 
ficients 



is called if dress is evalu- 
ated outside domain 



ansatz for dressing function 
with the transformation 
from cheb_f unc_traf 
taken into account, e.g., 
ghost dressing func- 
tion ansatz given by 
ansatzGh4dLog 
e.g., Ap^,p2 = Zi for pi,P2> 
(A/2)^ in scalar-gauge field 
vertex 



int 

*n_A 



numbers of expansion coeffi- 
cients in the dim_x variables 



numbers of expansion coeffi- 
cients for every external mo- 
mentum 



double 
*x_A 



(discrete) 
points 



interpolation 



external momenta (Matsub- 
ara frequencies) 



Member 



General purpose 



DSE usage/example 



Table C.7: User defined members of a quad structure. 
References to the scalar-gauge field vertex refer to the 
example of section [5} those to the ghost-gluon system to 
the example of section [6] 



Member 


General purpose 


DSE usage/example 


int 

bound_type 


defines if the integration 
boundaries depend on the 
nint_para sets of parame- 
ters 


defines if the integration 
boundaries depend on the 
external momenta 


Member 


General purpose 


DSE usage/example 
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Member 


General pm^pose 


DSE usage/example 


void 


defines the integration 


e.g., for the scalar-gauge 


(*boundary) 


boundaries for the different 


field vertex qsG [0, A^], ctl, 


(double *bound, 


integration regions 


ct2 e [-1, 1] 


double *x. 






int idim, void 






*int_parain) 






void 


constant factor of integrand 


constant factors of self- 


(*coeff ) 




energy kernels 


(double 






♦coefficients , 






void * int par am) 






int 


number of integration vari- 


number of loop momentum 


dim 


ables 


variables, e.g., 3 (qs, ctl 
and ct2) for scalar-gauge 
field vertex 


void 


initializes the nint_para 


automatically set to void 


(*init_para) 


parameter values 


outer .argument which ini- 


(int i, void 




tializes the external mo- 


*int_param) 




menta 


void 


the integrand 


kernels of the self-energies 


(* integrand) 






(double *erg, 






double *x, void 






*int_param) 






double 


Jacobian 


usually set to one 


(*jacob) 






(double *x) 






int 


number of different inte- 


number of dressing func- 


nint 


grands that are integrated 


tions multiplied by number 




over the same variables 


of graphs with the same in- 
tegration variables, e.g., 2 * 
2 = 4 for the scalar-gauge 
field vertex 


int 


number of times the inte- 


number of external mo- 


nint.para 


grands are integrated for 


menta at which the self- 




different parameter values 


energy is evaluated, i.e., 
number of interpolation 

points 


int 


number of integration re- 


e.g., the radial integration 


*n_part 


gions for any of the dim in- 


for the ghost-gluon system 




tegration variables 


is split into two parts 


int 


number of quadrature nodes 


e.g., the radial integration 


*nodes_part 


for the different integration 

regions 


of the ghost-gluon system 


Member 


General purpose 


DSE usage/example 
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Member 


General purpose 


DSE usage/example 


double 
*parain 


defines additional parame- 
ters in a quadrature rule, 
e.g., for the implemented 
double exponential 


not used in any of the 
DSE example, correct 
usage demonstrated in 
sphere^main. cpp 


int 

*traf 


defines the transformation 
function from [—1,1] to the 
true region of integration for 
different integration regions 


e.g., for the scalar-gauge 
field vertex the qs inte- 
gration is mapped via a 
modified logarithmic map- 
ping whereas the integra- 
tions in ctl and ct2 are not 
mapped 


int 
*type 


quadrature rules for the dif- 
ferent integration regions 


e.g., the scalar gluon vertex 
uses Fejer's second rule for 
the qs integration, Gauss- 
Chebyshev quadrature for 
the ctl integration and 
Gauss-Legendre quadrature 
for the ct2 integration 


Member 


General purpose 


DSE usage/example 



Table C.8: User defined members of a dse structure rele- 



vant for solving DSEs (not already included in table C.6). 



Member 


Short description 


double 


anomalous dimension of a dressing function; optional 


anoin_dini 




std: : string 


name of the DSE; optional 


DSE_name 




double 


stopping criterion for iteration: absolute difference of 


epsabs 


solutions 


double 


stopping criterion for iteration: relative difference of so- 


epsrel 


lutions 


double 


stepsize for Broyden's method; required only for secant 


h 


method 


void 


initializes coefficients *n_A and interpolation points, i.e.. 


(*init_dress) 


initial guess for the dressing function and definition of 


(void *dse_param) 


external momenta 


void 


initializes renormalization parameters 


(*init_renormparam) 




(void *dse_param) 




int 


counter for iterations; optional 


it.counter 




int 


stopping criterion for iteration: maximal iteration num- 


maxiter 


ber 


Member 


Short description 
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Member 


Short description 


void 


model parameters, usually the same for all Green func- 


(*niod) 


tions 


int 


number of diagrams with the same integral, *n_loop [i] , 


*n_loop 


iG {0, . . . ,n_loopnumber— 1} 


int 


number of different integrals 


n_loopnumber 




int 


number of other Green functions contributing to the 


n_otherGF 


given DSE 


struct dse 


one dse structure for every other Green function con- 


*otherGF 


tributing 


double 


stores the external momentum for a given index when 


* out er ^argument s 


the function outer_argument is called 


struct quad 


one quad structure for every of the n.loopnumber differ- 


*Q 


ent integrals 


void 


function that is called after calculating the loop integrals 


(*renorm) 


of the self-energy in iteration based solving routines 


(void *dse_param) 




int 


index of renormalization point 


*renorm_i 




int 


determines if *Z_renorin, *renorin_param, *renorm_x 


renorm.init 


and *renorm_i are allocated automatically by 
init_A_xA 


int 


number of renormalization parameters 


renormja_param 




int 


number of renormalization constants 


renornui_Z 




double 


stores renormalization parameters 


*renorin_parani 




double 


renormalization point 


*renorin_x 




double 


stores renormalization constants 


*Z_renorin 




Member 


Short description 
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